program afgeleide
    use, intrinsic IEEE_EXCEPTIONS
    implicit none
    real, parameter :: ALPHA = 3
    integer, parameter :: wp = kind(0e0)
    real(wp) :: startwaarde = 1.0
    real(wp) :: approx, exact
    integer :: i
    do i=1,16
        approx = approx_df(startwaarde,10.0_wp**(-i))
        exact = df(startwaarde)
        write(*, fmt='("h=10^",i3,", approx=",es9.2,", exact=",es9.2,"-> relerr=",es9.2)'),-i,approx,exact,(approx-exact)/exact
    enddo
contains

function f(x) result(value)
    real(wp), intent(in) :: x
    real(wp) :: value
    value = (x-ALPHA)**(-7)
end function

function approx_df(x,h) result(value)
    real(wp), intent(in) :: x, h
    real(wp) :: value
    value = (f(x+h)-f(x-h))/(2*h)
end function

function df(x) result(value)
    real(wp),intent(in) :: x
    real(wp) :: value
    value = -7*(x-ALPHA)**(-8)
end function

end program
    
